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INTRODUCTION: 


Palpation  by  clinicians  is  an  effective  examination  frequently  performed  for  breast  cancer  detection  and 
monitoring  treatment.  The  most  common  symptom  of  breast  cancer  is  a  lump,  and  studies  show  that  the 
majority  of  breast  cancers  are  found  by  palpation.  The  primary  goal  of  this  project  is  to  develop  a  tactile 
mapping  system  that  can  quantitatively  measure  the  location  and  applied  forces  in  breast  palpation,  and  the 
tactile  features  of  detected  breast  lumps.  By  combining  the  knowledge  of  vision  based  neural  networks,  tactile 
sensing  technologies  are  integrated  for  the  investigation  of  soft  tissue  interaction  with  tactile/force  sensors, 
where  the  hard  inclusion  (breast  lump)  can  be  characterized  through  neural  network  learning  capabilities, 
instead  of  using  a  simplified  complex  biomechanics  model  with  many  heuristic  assumptions.  Through  a 
nonlinear  soft  tissue  analysis  by  the  neural  network,  the  system  can  accurately  characterize  and  document  breast 
lumps.  These  maps  will  serve  as  an  objective  documentation  of  palpable  lesions  for  future  comparative 
examinations.  Limited  laboratory  experimental  evaluation  has  been  conducted  to  collect  preliminary  data  that 
will  form  a  full-scale  Task  II  research  plan.  Our  considerable  efforts  in  the  Task  II  are 

•  To  evaluate  die  sensitivity  and  reproducibility  of  the  TMD  system  based  on  the*  feedback  from  the 
model-based  laboratory  evaluation  and  preliminary  expert  evaluation. 

•  To  design  suitable  neural  network  based  intelligence  system  for  tactile  signal  interpretation 

•  To  develop  a  training  algorithms  and  computer  interpretation  codes. 
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BODY: 


In  the  first  year  of  this  award,  we  have  conducted  a 
systematic  study  to  adapt  new  tactile  sensing  technologies  to 
the  needs  of  improving  physical  breast  examination,  to  gather 
preliminary  data  regarding  potential  clinical  applications,  and 
to  advance  the  fundamental  understanding  of  palpation.  We 
established  that  an  existing  tactile  sensor  array  and  related 
force  measurement  technologies  could  be  modified 
specifically  for  use  as  a  tactile  mapping  device  (TMD)  and 
could  be  adapted  for  robotic  medical  sensing  applications 
[3,4,5].  We  applied  the  TMD  to  various  educational  models, 
which  allowed  the  collection  of  sufficient  tactile  images. 

In  the  second  year  of  this  award,  we  have  taken  the  feedback 
from  the  expert  evaluation  into  our  considerations  in 
designing  neural  network-based  intelligence  system.  Our  goal 
is  aimed  at  the  development  of  an  intelligence  system  to 
estimate  and  track  changes  in  tumors  over  time  during 
diagnosis  and  treatment.  The  precise  relationship  between 
tactile  quantities  and  the  breast  texture  and  abnormalities  is 
complex.  Our  approach  to  this  problem  is  to  build  on  tactile 
signal  processing  algorithms  developed  for  robotics 
applications,  which  relates  measured  surface  stresses  to 
shapes  of  contact.  To  analyze  tactile  abnormalities  with 
computational  intelligence,  we  have  designed  neural  network 
bases  to  accurately  estimate  and  track  changes  in  tumors  over 
time. 

The  precise  relationship  between  the  mounted  tactile 
quantities  and  the  breast  texture  and  abnormalities  is 
complex.  Some  work  has  recently  appeared  on  this  inversion 
problem  for  palpation  applications,  but  the  success  of  this 
approach  is  limited  because  it  relies  on  gross  force 
information  rather  than  the  pressure  distribution  information 
that  will  be  provided  by  our  tactile  array  sensors.  We  have 
previously  developed  signal  processing  algorithms  for  this 
purpose,  which  can  be  modified  to  map  the  sensor 
measurements  to  the  structure  of  the  breast  lump  [22,27]. 
This  approach  will  provide  a  better  understanding  of  breast 
palpation  and  establish  a  basis  for  constructing  refined 
palpation  devices.  Our  main  contribution  in  the  field  of 
medical  image  analysis  is  in  the  development  of  vision  type 
neural  networks,  fractal  imaging,  contextual  segmentation, 
and  advanced  pattern  recognition  techniques. 
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Figure  1.  Step  in  file  development  of 
TMD  system  for  breast  lump 
characterization  that  will  demonstrate  the 
feasibility  of  neural  network  classifies 
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We  have  put  considerable  effort  on  the  development  of  the  convolution  neural  network  (CNN)  and  its 
associated  training  methods.  Figure  3  shows  the  typical  structure  of  CNN.  The  structure  of  the  CNN  is  a 
simplified  version  of  “neocognitron,”  which  simulates  the  network  structure  and  signal  propagation  in 
vertebrate  animals.  The  structure  of  the  CNN  is  used  as  one  of  the  fundamental  research  tools  of  the  proposed 
research.  We  found  that  the  performance  of  the  CNN  is  more  powerful  than  a  regular  fully  connected  neural 
network,  such  as  the  Multilayer  perceptron  (MLP)  neural  network,  in  disease  pattern  recognition. 

The  Soft  tissue  analysis  for  the  characterization  of  palpable  breast  lesions  requires:  • 

(1)  The  quantitative  measurement  of  applied  forces  and  tactile  images  of  pressure  distribution  during 
palpation, 

(2)  The  construction  of  the  tactile  featured  database  and  training  algorithms,  and 

(3)  The  estimation  of  key  tactile  parameter  values  associated  with  detected  lesions  through  a  neural  network 
based  soft  tissue  modeling  and  analysis. 


InforratUte  3D  Rjunbcw  Cwaert 


Figure  2.  System  block  diagram  of  prototype  TMD. 

During  the  experiments,  the  forces/torques  applied  by  the  human  examiner  during  the  palpation  are  recorded  in 
terms  of  six  force  parameter  values:  Fx,  Fy,  Fz,  Mx,  My,  Mz.  The  TMD  is  used  to  acquire  various  tactile 
images,  by  palpating  the  site  after  the  examiner  locates  the  lesion  by  initial  hand  palpation.  Since  the 
relationships  between  the  tactile  images  and  lesion  characteristics  are  expected  to  be  complex  and  nonlinear,  we 
believe  that  the  inverse  problem  (extraction  of  lesion  characteristics  from  tactile  images)  can  be  solved  when 
sufficient  tactile  information  is  acquired.  As  a  result,  for  each  lesion,  after  pushing  the  TMD  against  the  breast 
to  achieve  certain  level  of  force,  the  examiner  rotates  the  TMD  to  five  different  orientations  at  each  of  which  the 
information  about  the  forces/torques  and  tactile  image  are  simultaneously  recorded.  OR1  is  the  initial 
orientation  of  the  TMD  after  being  pushed  by  the  examiner  to  a  certain  level  of  force.  OR2,  OR3,  OR4  and  OR5 
are  the  TMD  orientations  after  it  is  rotated  forward,  backward,  left  and  right,  respectively.  Figure  2.  shows  the 
prototype  of  the  tactile  mapping  system,  where  an  8x8  tactile  sensor  array  and  a  6DOF  force/torque  sensor  is 
employed.  In  order  to  demonstrate  the  effectiveness  of  the  tactile  mapping  device  for  the  characterization  and 
documentation  of  detected  lumps,  intensive  experiments  have  been  conducted  to  pilot-test  the  sensitivity  and 
reproducibility  of  the  system  in  measuring  the  applied  forces  and  pressure  profiles  due  to  lumps.  Examining  the 
graphical  displays  in  Figure  5.b  and  5.c  in  comparison  with  that  of  the  case  with  no  lesions  (Figure  5.a),  we 
noticed  that  a  peak  consistently  occurred  in  the  display.  We  can  draw  preliminary  conclusions  that  a  peak  in  the 
graphical  display  indicates  the  existence  of  a  lesion.  The  output  from  our  TMD  will  be  used  as  an  input  for  our 
neural  network. 
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A  convolution  neural  network  (CNN)  based  soft  tissue  modeling  and  analysis  technique  is  developed  to 
estimate  lesion  sizes  and  depth,  shown  in  Figure  4,  where  each  input  denotes  a  cluster  corresponding  to  the 
tactile  image  of  one  force  projection,  and  each  hidden  cluster  extracts  appropriate  tactile  features  to  highlight 
mass  characteristics.  The  CNN  outputs  consist  of  two  sets  of  signals:  parameter  estimates,  and  reconstructed 
tactile  profiles. 

In  neural  network  based  mass  detection,  we  collected  the  output  from  our  TMD  used  the  output  as  inputs  to  our 
neural  network.  We  separated  our  database  into  two  sets;  one  set  will  be  used  for  training,  and  the  remaining  set 
will  be  used  for  testing.  This  database  was  used  to  train  two  neural  network  systems: 

(1)  A  conventional  3 -layer  back  propagation  neural  network 

(2)  A  CNN  training  method 

The  error  back  propagation  algorithm  trained  both  neural  network  systems  by  feeding  the  features  from  the 
input  layer,  and  registering  the  corresponding  target  value  at  the  output  side.  Once  the  training  of  the  neural 
networks  is  complete,  we  will  then  use  the  remaining  sets  for  the  testing. 


.  .  *  .Tavvrju  j 

Cluster  Output  IMAGE 


Figure  3.  Convolution  neural  network  (CNN)  based  breast  tumor  characterization  and  parameter  estimation 
where  the  inputs  are  the  tactile  images  and  associated  applied  angles  and  forces,  and  the  outputs  are  the 
estimated  size  and  depth  of  the  lump  and  associated  probability  of  being  a  tumor. 
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Figure  4.  Shows  the  outputs  from  our  TMD  will  be  used  as  an  input  to  our  convolution  neural  network  based 
breast  tumor  characterization  and  parameter  estimation. 
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We  applied  several  image  processing  and  training  methods  designed  to  mimic  a  physician’s  interpretation  using 
the  breast  lump  database.  The  TMD  will  be  evaluated  by  physicians,  and  will  undergo  a  limited  evaluation  of 
diagnosis  and  documentation  using  simulations.  Clinical  evaluations  will  be  carried  out  under  the  supervision 
of  Dr.  Matthew  Freedman.  In  particular,  we  will  evaluate  the  sensitivity,  accuracy,  and  precision  of  the  tactile 
sensors  during  simulated  CBE  applications. 


KEY  RESEARCH  ACCOMPLISHMENTS: 

Our  key  goal  of  Task  II  is  to  develop  a  neural  network  based  intelligence  system  that  estimates  and  tracks  the 
changes  (in  size  and  depth)  of  the  tumors  over  time  after  diagnosis  and  during  treatment  to  develop  a  neural 
network.  Our  accomplishments  in  Task  n  are  the  following: 

•  By  taking  the  feedback  from  experts,  we  have  evaluated  the  TMD’s  performance  in  improving  diagnosis 
in  terms  of  precision,  sensitivity,  and  specificity.  We  have  re-designed  a  prototype  TMD  system  based 
on  the  feedback  from  the  model-based  laboratory  evaluations,  and  preliminary  expert  evaluations 

•  We  have  designed  a  suitable  neural  network-based  intelligence  system  for  tactile  signal  interpretation. 
By  combining  the  knowledge  of  vision-based  neural  networks,  tactile  sensing  technologies  are 
integrated  for  the  investigation  of  soft  tissue  interactions  with  the  tactile/force  sensor,  where  hard 
inclusions  (breast  lump)  can  be  characterized  through  neural  network  learning  capabilities 

•  We  have  developed  neural  network  training  algorithms  and  computer  interpretation  codes,  convolution 
neural  network  (CNN)  based  breast  tumor  characterization  and  parameter  estimation;  the  inputs  are  the 
tactile  images  and  associated  applied  angles  and  forces;  the  outputs  are  the  estimated  size  and  depth  of 
the  lump  and  associated  probability  of  being  a  tumor. 
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REPORTABLE  OUTCOMES: 

•  R.  Srikanchana,  K.  Huang,  J.  Xuan,  M.  Freedman,  and  Y.  Wang,  “Mixture  of  Principal  Axes 
Registration  for  Change  Analysis  in  Computer-Aided  Diagnosis,”  30th  Applied  Imagery  Pattern 
Recognition  Workshop,  Washington,  DC,  October  2001. 

•  R.  Srikanchana,  Y.  Wang,  C.  C.  Nguyen,  and  M.  Freedman,  “Tactile  Imaging  of  Palpable  Breast 
Cancer,”  Proc.  SPIE  Medical  Imaging,  vol.  4682,  San  Diego,  CA,  February  2002. 

•  R.  Srikanchana,  J.  Xuan,  M.  Freedman,  and  Y.  Wang  “Mixture  of  principal  axes  registration:  a  neural 
computation  approach,”  Proc.  SPIE  Medical  Imaging,  vol.  4684,  San  Diego,  CA,  February  2002 

•  Presented  the  paper  Mixture  of  Principal  Axes  Registration  for  Change  Analysis  in  Computer-Aided 
Diagnosis,”  3Cf  Applied  Imagery  Pattern  Recognition  Workshop ,  Washington,  DC,  October  2001 . 

•  Presented  the  poster  “Tactile  Mapping  and  Tumor  Modeling  for  Breast  Cancer  Detection,”  The  third 
Era  of  Hope  meeting  for  the  Department  of  Defense  (DOD)  Breast  Cancer  Research  Program  (BCRP ) 
Orlando,  FL,  September  2002. 


CONCLUSIONS: 

This  project  is  concerned  with  the  detection  of  lumps  in  breast  palpation  using  a  tactile  mapping  device  (TMD) 
system.  It  presented  the  main  components  of  the  TMD  and  explained  their  operations.  Experiments  were 
conducted  on  a  breast  model  to  detect  different  types  of  lumps.  Experimental  data  showed  that  the  peak  in  the 
display  of  pressure  distribution  of  the  TMD  indicated  the  presence  of  a  lump.  The  knowledge  of  vision-based 
neural  networks  and  tactile  sensing  technology  are  integrated  for  the  first  time  for  investigations  of  soft  tissue 
interactions  with  the  tactile/force  sensor;  the  hard  inclusion  (breast  cancer)  can  be  characterized  through  neural 
network  learning  capabilities,  instead  of  using  a  simplified  complex  biomechanics  model  with  many  heuristic 
assumptions.  Tactile  mapping  systems  using  convolution  neural  networks  and  tactile  sensing  arrays  can  extract 
invariant  properties  of  the  detected  breast  lump,  and  make  it  possible  for  the  first  time  to  quantitatively  and 
objectively  record  the  processes  and  findings  of  a  breast  palpation.  With  the  proven  power  of  nonlinear  signal 
processing  using  both  convolution  neural  networks  (CNN)  and  multi-layer  perceptron  (MLP),  we  expect  that 
the  tactile  parameters  associated  with  the  lesions  (i.e.,  the  size  and  depth  of  the  lesion)  can  be  estimated  more 
accurately  than  those  by  the  conventional  approaches  using  the  first  principle  of  engineering. 
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Abstract 


Non-rigid  image  registration  is  a  prerequisite  for  many 
medical  image  analysis  applications  such  as  image  fusion 
of  multi-modality  images  and  quantitative  change  analysis 
of  a  temporal  sequence  in  computer-aided  diagnosis .  By  es¬ 
tablishing  the  point  correspondence  of  the  extracted  feature 
points t  it  is  possible  to  recover  the  deformation  using  non¬ 
linear  interpolation  methods  such  as  die  thin-plate-spline 
approach  However  it  is  a  difficulty  task  to  establish  exact 
point  correspondence  due  to  the  high  complexity  of  nonlin¬ 
ear  deformation  existed  in  medical  images.  In  this  paper, 
a  mixture  of  principal  axes  registration  (mPAR)  method  is 
proposed  to  resolve  the  correspondence  problem  through  a 
neural  computational  approach.  The  novel  feature  ofmPAR 
is  to  align  two  point  sets  without  needing  to  establish  ex¬ 
plicit  point  correspondence .  Instead ,  it  aligns  the  two  point 
sets  by  minimizing  the  relative  entropy  between  their  prob¬ 
ability  distributions  resulting  in  a  maximum  likelihood  es¬ 
timate  of  the  transformation  matrix.  The  registration  pro¬ 
cess  consists  of:  (1)  a  finite  mixture  scheme  to  establish  an 
improved  point  correspondence  and  (2)  a  multilayer  per- 
ceptron  neural  network  (MLP)  to  recover  the  nonlinear  de¬ 
formation.  The  neural  computation  for  registration  used  a 
committee  machine  to  obtain  a  mixture  of  piece-wise  rigid 
registrations ,  which  gives  a  reliable  point  correspondence 
using  multiple  extracted  objects  in  a  finite  mixture  scheme. 
Then  the  MLP  is  used  to  determine  the  coefficients  of a  poly¬ 
nomial  transform  using  extracted  cross  points  of  elongated 
structures  as  control  points.  We  have  applied  our  mPAR 
method  to  a  temporal  sequence  of  mammograms  of  a  single 
patient.  The  experimental  results  show  that  mPAR  not  only 
improves  the  accuracy  of  the  point  correspondence  but  also 
results  in  a  desirable  error-resilience  property  for  control 
point  selection  errors. 


1  Introduction 

Image  registration  is  an  essential  step  for  many  medical 
image  analysis  applications  such  as  image  fusion  of  multi- 
modality  images  and  quantitative  change  analysis  of  a  tem¬ 
poral  sequence  in  computer-aided  diagnosis.  Due  to  com¬ 
plicated  by  misregistrations  caused  by  patient  motion  and 
physical  change,  the  estimation  of  transformational  geom¬ 
etry  from  two  point  sets  is  an  essential  step  to  medical  im¬ 
age  fusion  and  computer  vision  [1,  2].  The  task  is  to  re¬ 
cover  a  matrix  representation  requiring  a  set  of  correspon¬ 
dence  matches  between  features  in  the  two  coordinate  sys¬ 
tem  [3].  Assume  two  3-D  data  point  sets  {p^}  and  {piB  }; 
i  =  1?  2,  N  are  related  by 

PiB  —  R PiA  +  T  -I-  N»  (1) 

where  R  is  a  rotation  matrix,  T  is  a  translation  vector,  and 
N i  is  a  noise  vector*  Given  {piA}  and  {pis},  Arun  et.  al 
present  an  algorithm  for  finding  die  least-squares  solution 
of  R  and  T,  which  is  based  on  the  decoupling  translation 
and  rotation  and  the  singular  value  decomposition  (SVD)  of 
a  3  x  3  cross-covariance  matrix  [3]* 

The  major  limitation  of  the  present  method  is  twofold: 

(1)  while  feature  matching  methods  can  give  quite  accu¬ 
rate  solutions,  obtaining  correct  correspondences  of  fea¬ 
tures  is  a  hard  problem,  especially  in  the  cases  of  images 
acquired  using  different  modalities  or  from  inter-subjects; 
and  (2)  a  rigidity  assumption  is  heuristically  imposed  while 
those  newly  developed  deformable  matching  methods  can 
be  computationally  expensive  and  typically  need  good  ini¬ 
tial  guesses  to  assure  correct  convergence  [2, 4],  One  pop¬ 
ular  method  that  is  correspondenceless  is  the  principal  axes 
registration  (PAR)  [1],  that  is  based  on  the  relatively  stable 
geometric  properties  of  image  features,  ie,,  the  geometric 
information  contained  in  these  stable  image  features  is  often 
sufficient  to  determine  the  transformation  between  images 

[2] ,  Once  again,  it  cannot  handle  the  cases  with  non-rigid 
objects* 

We  present  a  neural  computation  based  non-rigid  reg¬ 
istration  using  piece-wise  rigid  transformation*  The  novel 


1 


feature  is  to  align  two  point  sets  without  needing  to  estab¬ 
lish  explicit  point  correspondences,  where  the  derivation  is 
realized  by  minimizing  the  relative  entropy  between  the  two 
point  distributions  resulting  in  a  maximum  likelihood  esti¬ 
mate  of  the  transformation  matrix.  That  is  rather  than  using 
a  single  transformation  matrix  which  give  a  large  registra¬ 
tion  error,  we  attempt  to  inteipolatively  apply  a  mixture  of 
transformations.  By  further  generalizing  PAR  to  a  mixture 
of  principal  axes  registration  (mPAR)  scheme,  with  soft  par¬ 
titioning  of  the  data  set,  the  mixture  is  fit  using  expectation- 
maximization  algorithm,  to  estimate  the  transformational  of 
the  orthogonal  set  of  eigenvalues  and  eigenvectors  of  the 
auto-covariance  matrix.  By  applying  a  committee  machine 
to  a  non-rigid  registration,  we  can  acquired  the  registration 
based  on  a  mixture  of  piece-wise  transformations  of  the  data 
set.  Then  the  correspondences  control  points  are  obtained. 
As  a  final  step,  the  warped  image  is  obtaining  using  the  neu¬ 
ral  network  based  non-linear  mapping,  to  obtain  the  poly¬ 
nomial  transform  based  on  extracted  control  points  using 
multilayer  perception  neural  network. 

2  Method 

Suggested  by  information  theory  [5],  we  note  that,  since 
the  control  point  sets  in  two  images  can  be  considered  as 
two  separate  realizations  of  the  same  random  source,  we 
do  not  need  to  establish  point  correspondences  to  extract 
the  transformation  matrix.  In  other  words,  if  we  denote  by 
P{Pi  }  the  distribution  of  the  control  point  set  in  an  image, 
we  have  the  simple  relationship 

P{pjb}  =  ^{Rp^-fT}  +  €  (2) 

•Jk  ' 

where  e  is  the  noise  component  Since  the  probability  dis¬ 
tributions  can  be  computed  independently  on  each  image 
without  any  need  to  establish  feature  correspondences,  and 
given  the  two  distributions  of  the  control  point  sets  in  the 
two  unages,  we  can  recover  the  transformation  matrix  in  a 
simple  fashion  [2]. In  order  to  establish  feature  correspon¬ 
dence  between  the  two  images,  we  applied  our  segmenta¬ 
tion  method  to  extract  the  control  points  between  two  im¬ 
ages.  After  the  control  has  been  established,  we  applied  our 
registration  method  to  the  images,  as  we  now  sketch. 

2.1  Control  Point  Extraction 

Assume  that  each  pixel  in  the  MR  image  can  be  decom¬ 
posed  into  pixel  images  x  and  context  image  L  By  ignoring 
information  regarding  the  spatial  ordering  of  pixels,  we  can 
treat  context  images  as  random  variables  and  describe  them 
using  a  multinomial  distribution  with  unknown  parameters 
tt/c,  &  =  1, AT,  which  can  be  interpreted  as  a  prior  proba¬ 
bility  of  pixel  labels  determined  by  the  global  context  infor¬ 
mation.  In  particular,  based  on  the  statistical  properties  of 


MR  pixel  images,  where  pixel  image  is  defined  as  the  ob¬ 
served  gray  level  associated  with  the  pixel,  use  of  Standard 
Finite  Normal  Mixture  (SFNM)  distribution  is  justified  to 
model  the  image  histogram  by  determining  the  optimal  pa¬ 
rameters  with  respect  to  a  distance  measure  of  a  sum  of  the 
following  general  form: 


K 

Mx)  =  '*Tirkg(x\fik,ol) 


k=l 


K 

with  nk  =  1, 7T*  >  0,  and 
i 


where  pk  and  ak  are  the  mean  and  variance  of  the  jfcth 
Gaussian  kernel,  and  7rk  is  the  global  regularization  pa¬ 
rameter.  We  use  K  to  denote  the  number  of  Gaussian 
components  and  r  to  denote  the  parameter  vector.  This 
tissue  quantification  is  achieved  through  three  completely 
unsupervised  steps:  (1)  parameter  initialization  by  optimal 
histogram  quantization;  (2)  model  estimation  by  histogram 
based  fast  expectation-maximization  (EM)  algorithm;  and 
(3)  model  selection  by  minimum  description  length  (MDL) 
criterion.  It  is  shown  that  the  SFNM  model  converges  to 
the  true  distribution  when  the  pixel  image  are  asymptoti¬ 
cally  independent  [18]. 

After  we  obtain  the  optimal  parameters  for  all  compo¬ 
nents,  we  use  a  multiple  thresholding  procedure  for  ini¬ 
tializing  image  segmentation  based  on  maximum  likelihood 
(ML)  principle,  which  is  followed  by  contextual  Bayes  re¬ 
laxation  labeling  (CBRL)  algorithm  to  obtain  a  consistent 
labeling  solution  based  on  localized  SFNM  formulation  for 
moproving  initial  segmentation  by  using  neighborhood  text 
regularities. 

We  define  the  component  in  localized  SFNM  by  the  sup¬ 
port  function: 


S<(k)  =  *® 


%/2 ITCTk 


exp(- 


) 


where  tt^  is  the  local  conditional  prior  of  regions,  the 
support  function  Si(k)  is  a  function  of  the  component  (tis¬ 
sue  type)  k .  Then  tissue  segmentation  is  interpreted  as  the 
satisfaction  of  a  system  of  inequalities  as  follows: 


Si(k)  >  Si(k ) 

for  all  k  and  for  i  =  1,...,#,  where  liyi  =  1, N  is  the 
context  images,  and  a  consistent  labeling  is  defined  as  the 
one  having  maximum  support  at  each  pixel  simultaneously. 
We  further  define  the  average  local  consistency  measure: 


N 

i=l  k 

to  Me  consistent  labeling  to  global  optiinization[19]. 

2.2  Registration 

For  observation  of  tbe  distributions,  we  can  estimate  R 

distance)  between  P{PjB)  and  P{rPu+T}.  Tbe  least  rela¬ 
tive  entropy  estimator  is  then  defined  as 

I  IP{Rp*+T})  (3) 

where  D  denotes  the  relative  entropy  measure.  Following 
the  same  strategy  to  decoupling  translation  and  rotation  as 
in  [3],  we  can  define  a  new  data  point  by  —  piA  - 
and  qjs  =  Pjb  —p% ,  where  p°A  and  are  the  centroids  of 
{PiA}  and  {pj£ }  respectively.  Then  the  optimal  geometric 
transformations,  R  and  T,  are  defined  as 

R  =  UbHU^  andT  =  p|-Rp*  (4) 

where  the  superscript  t  denotes  matrix  transposition,  XJA 
and  UB  are  3  x  3  orthonormal  matrices,  and  H  is  a  3x3 
diagonal  matrix  with  element  hm  =  ^/XmB/XmA.  Note 
that  the  transformation  U  consists  of  the  orthonormal  set  of 
eigenvectors  and  hm  is  the  squared  root  of  die  eigenvalues 
Am  of  the  auto-covariance  matrix  C  for  m  =  1, 2, 3  and  for 
{piA}  and  {pjB },  respectively. 

However,  because  of  its  global  linearity,  the  application 
of  PAR  is  necessarily  somewhat  limited  [6],  An  alternative 
paradigm  is  to  model  a  multimodal  control  point  set  with  a 
collection  of  local  linear  models.  The  method  is  a  two-stage 
procedure:  a  soft  partitioning  of  the  data  set  followed  by  es¬ 
timation  of  the  principal  axes  within  each  partition  [7].  Re¬ 
cently  there  has  been  considerable  success  in  using  standard 
finite  normal  mixture  (SFNM)  to  model  die  distribution  of 
a  multimodal  data  set  and  the  association  of  a  SFNM  dis¬ 
tribution  with  PAR  offers  the  possibility  of  being  able  to 
register  two  images  through  a  mixture  of  probabilistic  prin¬ 
cipal  axes  transformations  [7].  Assume  that  there  are  Kq 
control  point  clusters,  where  each  control  point  cluster  de¬ 
fines  a  transformation  {R*5  T&}.  Thus  for  a  point  p^,  its 
new  locations  corresponding  to  each  of  the  transformations 
are  p nk  =  R kPnA  4-  Tk  for  k  =  1,  Further  assume 

that  the  control  point  set  defines  a  SFNM  distribution 

Xo 

/(Pi)  =  ]TV  *£(PilMk, C*)  <5> 

k=i 

where  g  is  the  Gaussian  kernel  with  mean  vector  and 
auto-covariance  matrix  C&,  and  wk  is  the  mixing  factor  pro¬ 
portional  to  the  number  of  control  points  in  cluster  k.  For 


each  of  the  control  point  sets  {piA  }  and  {piB  }?  the  mixture 
is  fit  using  the  expectation-maximization  (EM)  algorithm. 
The  E  step  involves  assigning  to  tire  linear  models  respon¬ 
sibilities  from  the  control  points;  the  M  step  involves  re- 
estimating  the  parameters  of  the  linear  models  in  die  light 
of  this  assignment  [7]. 

Thus  the  statistical  membership  of  point  pnA  belonging 
to  each  of  the  control  point  clusters  can  be  derived  by 

■Snfc  =  P(Rk,Tk\j>nA)  =  ^A9(pnA^kA-CkA) 

f{PnA ) 

i.e.,  the  posterior  probability  of  {R*,  Tfc}  given  pnA.  Thus, 
we  can  define  the  FMR  transformation  as 


p»=E  = £  ^<py'^'c“)(RtPnJ+Tt) 


where  {Rk,Tk}  is  determined  based  on  f(pi)  = 
12k=i  KkdiPilP'ki  Cjfe)  that  we  have  estimated  in  the  pre¬ 
vious  step  using  the  EM  algorithm.  Note  that  now  we  do 


need  die  correspondences  between  die  two  control  point 
clusters,  and  these  correspondences  may  be  found,  after  a 


global  PAR  is  initially  performed,  by  using  a  site  model 


supported  approach,  or  a  dual-step  EM  algorithm  to  unify 
the  teste  of  estimating  transformation  geometry  and  identi¬ 


fying  cluster-correspondence  matches  [4].  This  philosophy 
for  recovering  transformational  geometry  of  the  non-rigid 
objects  is  similar  in  spirit  to  the  divide-and-conquer  princi¬ 
ple  [6],  under  which  die  relative  entropy  between  the  two 


point  sets  reaches  its  minimum 


arES.z,(^'*>l|p<ES1««p.p.+T.)>) 


both  globally  and  locally. 

Based  on  a  mixture  of  probabilistic  principal  axes  trans¬ 
formations,  the  next  section  describes  a  neural  computation 
using  a  committee  machine  approach  for  which  a  complex 
computational  task  is  solved  by  dividing  it  into  a  number  of 
computationally  simple  taste  and  then  combining  the  solu¬ 
tions  to  those  taste. 


2.2.1  Neural  Computation 

A  neural  network  interpretation  of  the  EM  algorithm  is 
given  in  [8].  Because  of  its  reputation  of  being  slow  in 
which  new  information  acquired  in  the  expectation  step  is 
not  used  immediately,  on-line  versions  of  the  EM  algorithm 
are  proposed  for  large-scale  sequential  learning.  Thus,  we 
adopt  a  fully  unsupervised  and  incremental  stochastic  learn¬ 
ing  algorithm.  The  scheme  provides  winner-takes-in  prob¬ 
ability  (Bayesian  “soft”)  splits  of  the  control  points,  hence 
allowing  the  data  to  contribute  simulteneously  to  multiple 
clusters  which  results  in 


2(i+l)fc  = 

M  -  Step 


^JL°g(p*+xi#*£°,  c£*>) 


(9) 


4<+1)-“#*i° +«(<)(!>*+! 


(10) 


Ci<+1)  =  ci<)+6(i)[(pi+a-MW)(pi+1-^))r_c(Ok<+i)fc 

^+1)  =  Th^)+jhz^k  (12) 


for  A;  =  and  for  {p^}  and  {pis},  respectively, 

where  a(i)  and  b(i )  are  introduced  as  the  learning  rates,  two 
sequences  converging  to  zero,  ensuring  unbiased  estimates 
after  convergence.  This  procedure  is  termed  as  neural  com¬ 
putation  of  the  EM  algorithm,  where  at  each  complete  cycle 
of  the  algorithm,  we  first  use  “old!'  set  of  parameter  val¬ 
ues  to  determine  the  posterior  probabilities  Z(<+ i)fc.  These 
posterior  probabilities  are  then  used  to  obtain  “new"  val- 
ues  tt£+: 1\  /a£+1\  The  algorithm  cycles  back  and 

forth  until  the  value  of  relative  entropy  between  the  data 
histogram  and  mixture  model  reaches  its  minimum 


Figure  1.  A  pair  of  real  mammograms  taken 
over  a  period  of  time  from  the  same  patient. 


are  obtained.  The  final  step  is  to  calculate  the  polynomial 
transformation  using  piece-wise  interpolation.  In  our  case 
for  a  given  corresponding  control  points  pair  of  the  image, 
we  considered  neural  network  based  MLP  to  recover  the 
nonlinear  deformation.  The  final  result  shows  the  warped 
image  that  has  been  corrected  by  applied  our  neural  com¬ 
putation  to  correct  most  of  the  scale  different  between  the 
images. 


Tc  ^{pjll/tPi))  (13) 

for  {piA}  and  {p^},  respectively. 

With  a  soft  partitioning  of  the  data  set  using  Eqs.  (9-12), 
control  points  will  now  effectively  belong  to  more  tha-n  one 
cluster  spatially.  Thus,  the  effective  input  values  are  p#  = 
*ik  (p i  —f^k)  an  independent  registration  transformation 
fc  in  the  committee  machine  [9].  We  then  extend  our  adap¬ 
tive  principal  components  extraction  (APEX)  algorithm  to 
a  probabilistic  version,  i.e.,  PAPEX  [10],  to  determine  Ufc 
for  {ptx}  and  {pis}>  respectively. 

Consider  the  modular  network  nature  of  Eq.  (7),  it  is  a 
mixture  of  local  expert  model,  in  which  the  individual  re¬ 
sponses  of  the  experts  are  nonlinearly  combined  by  means 
of  a  single  gating  network.  It  is  assume  here  that  the  differ¬ 
ent  experts  work  best  in  different  regions  of  the  input  space 
in  accordance  with  the  probabilistic  generative  model.  In 
our  case,  the  effective  input  values  are  p**  =  zik  (p*  -  fik) 
for  an  independent  registration  transformation  k  in  the  net¬ 
works.  A  committee  machine  consists  of  k  supervised  mod¬ 
ules,  with  a  soft  partitioning  of  the  data  set  using  EM  al¬ 
gorithm  called  experts,  and  an  integrating  unit  called  a  gat¬ 
ing  network  that  performs  the  function  PAPEX  to  determine 
the  transformation  of  the  orthogonal  set  of  eigenvalues  and 
eigenvectors  of  the  auto-covariance  matrix  among  the  ex¬ 
pert  networks.  The  output  of  our  committee  machine  is  a 
transformational  matrix  of  image  pair.  Based  on  those  cor¬ 
respondences  between  the  two  images,  the  control  points 


3  Results  and  Discussion 

We  applied  our  algorithm  to  a  temporal  sequence  of 
mammograms  of  a  single  patient  [15].  In  this  example, 
the  committee  machine  is  used  to  obtain  an  initial  registra¬ 
tion  using  multiple  extracted  objects  (skinline,  dense  tissue 
regions)  in  a  finite  mixture  scheme.  Then  MLP  was  used  to 
determine  the  coefficients  of  a  polynomial  transform  using 
extracted  vertical  and  horizontal  cross  points  of  elongated 
structures  as  control  points.  Previously,  thin-plate  spline 
(TPS)  was  used  to  determine  the  transform  coefficients  [13]. 
As  a  comparison,  we  consider  both  results  here.  Figure  1 
shows  the  raw  sequence  and  Figure  2  shows  the  resulting 
warped  image  for  MLP  (left)  and  TPS  (right).  Using  im¬ 
age  difference  (Figure  3)  we  see  that  the  most  of  the  scale 
different  between  images  has  been  corrected  by  using  MLP. 
While  the  TPS  distorts  the  image,  this  distortion  can  be  at¬ 
tributed  to  control  point  selection  and  correspondence.  The 
MLP  better  adapts  to  the  error  present  in  the  control  points 
thus  yielding  a  smoother  result. 

We  applied  our  approach  for  the  non-rigid  registration 
of  contrast-enhanced  breast  MRI.  A  3D  breast  MRI  scan  is 
acquired  prior  to  the  injection  of  contrast  agent,  followed 
by  a  sequence  of  3D  breast  MRI  scans  after  the  agent  has 
been  applied.  The  goal  of  image  registration  in  contrast- 
enhanced  breast  MRI  is  to  relate  any  point  in  the  post¬ 
contrast  enhanced  sequence  to  the  pre-contrast  enhanced 


Figure  2.  The  result  of  warped  current  image 
to  previous  image  (left:  MLP;  right:  TPS). 


Figure  4.  Segmentation  of  Fibre-Tissue 


Figure  5.  A  pair  of  a  contrast-enhanced  in 
MR!  (Left:  pre-contrast  image;  right:  post¬ 
contrast  image.) 


Figure  3.  Image  difference  (left:  MLP;  right: 
TPS) 


reference  image.  After  we  performed  a  pre-registration  us¬ 
ing  PAR  method  to  correct  the  global  alignment  between 
images,  we  applied  our  segmentation  algorithm  to  extract 
the  control  objects  between  two  images.  Figure  4  shows 
the  segmentation  of  Fibre-tissue.  Next,  we  applied  the  pro¬ 
posed  non-rigid  registration  to  these  images  using  the  con¬ 
trol  objects  that  has  been  acquired  by  the  segmentation.  Fig¬ 
ure  5  shows  a  pre-  and  post-contrast  enhanced  image  of  the 
patient  data  set  with  a  propose  registration  method.  Figure 
6  shows  the  post-contrast  enhanced  image  and  correspond¬ 
ing  difference  images  after  registration.  The  results  demon¬ 
strate  that  the  simple  image  differencing  can  isolate  such 
area. 

4  Conclusions 

In  this  paper,  we  presented  the  theoretical  concepts  and 
methods  of  a  neural  computation  based  non-rigid  registra¬ 
tion  algorithm  for  computer-aided  diagnosis.  The  approach 
uses  a  committee  machine  to  recover  the  total  transforma¬ 
tional  geometry  of  the  non-rigid  object  using  multiple  rigid 
transforms  combined  together  in  a  mixture  of  principal  axes 
registration  scheme.  Finite  mixture  transform  combination 
is  a  novel  technique  to  combine  multiple  transforms  that  are 
contained  in  a  single  image.  Other  than  local  transforms  no 
other  method  combines  multiple  transforms  together.  In  ad- 


Figure  6.  Difference  image  to  highlight  en¬ 
hanced  area. 


dition  finite  mixture  combinations,  yields  a  smooth  image 
while  local  transforms  yield  an  image  containing  discon¬ 
tinues  on  transform  boundaries.  In  addition,  the  registra¬ 
tion  obtained  in  the  committee  machine  is  fine  tuned  using 
a  non-linear  transform  generated  by  a  MLP  network  using 
extracted  control  points. 

We  applied  our  non-rigid  algorithm  to  a  temporal  se¬ 
quence  of  breast  image  of  a  single  patient  Some  distortion 
can  be  seen  in  the  final  warped  images  because  of  the  er¬ 
ror  in  control  point  selection  and  correspondence.  Improve¬ 
ment  is  this  portion  should  decrease  distortion  and  yield 
a  smoother  looking  image.  Using  neural  networks  in  this 
problem  has  increased  the  generality  of  this  approach  by 
allowing  the  algorithm  to  adjust  performance  as  imaging 
condition  change. 
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ABSTRACT 

This  paper  presents  the  development  of  a  prototype  Tactile  Mapping  Device  (TMD)  system  comprised  mainly 
of  a  tactile  sensor  array  probe  (TSAP),  a  3-D  camera,  and  a  foree/torque  sensor,  which  can  provide  the  means 
to  produce  tactile  maps  of  the  breast  lumps  during  a  breast  palpation.  Focusing  on  the  key  tactile  topology 
features  for  breast  palpation  such  as  spatial  location,  size/shape  of  the  detected  lesion,  and  the  force  levels  used 
to  demonstrate  the  palpable  abnormalities,  these  maps  can  record  the  results  of  clinical  breast  examination 
with  a  set  of  pressure  distribution  profiles  and  force  sensor  measurements  due  to  detected  lesion.  By  combining 
the  knowledge  of  vision  based,  neural  networks  and  tactile  sensing  technology;  the  TMD  is  integrated  for  the 
investigation  of  soft  tissue  interaction  with  tactile/force  sensor,  where  the  hard  inclusion  (breast  cancer)  can 
be  characterized  through  neural  network  learning  capability,  instead  of  using  simplified  complex  ^m^nice 
model  with  many  heuristic  assumptions.  These  maps  will  serve  as  an  objective  documentation  of  palpable  lesions 
for  future  comparative  examinations.  Preliminary  results  of  simulated  experiments  and  limited  pre-dinica! 
evaluations  of  the  TMD  prototype  have  tested  this  hypothesis  and  provided  solid  promising  data  showing  the 
feasibility  of  the  TMD  in  real  clinical  applications. 

Keywords:  T  actile  Mapping  Device,  Tactile  Sensor  Array  Probe,  force/torque  sensor,  breast  palpation,  breast 
cancer,  tactile  imaging. 


1.  INTRODUCTION 

Physical  breast  examination  is  an  effective  and  completely  non-invasive  method  for  the  detection  of  breast 
cancer.  *  With  a  lump  as  the  most  common  symptom  of  breast  cancer,  studies  show  that  the  majority  of  breast 
cancers  was  found  by  palpation  which  complements  mammography,  since  palpation  can  evaluate  breast  tissue 
near  the  chest  wall  and  axilia  that  is  not  accessible  to  mammography.3  In  addition,  studies  have  found  that  as 
many  as.  12-15%  of  cancers  that  were  detected  by  physical  examination  were  not  apparent  on  mammograms.1*3 

Unfortunately,  breast  palpation  has  been  hampered  by  problems  inherent  in  its  subjective  nature,  leading 
to  difficulty  in  interpreting  and  documenting  the  examiner’s  impressions  of  the  perceived  lump  in  terms  of 
tumor  characteristics.3  For  example,  a  physician  may  determine  that  a  palpable  suspicious  abnormality  needs 
continued  monitoring.  This  requires  maintaining  a  record  of  the  examination  results,  which  at  present  is  limited 
to  verbal  notes  about  parameters  such  as  the  position,  size,  and  hardness  of  the  lump.  Because  it  is  difficult  to 
verbalize  tactile  sensations,  the  subjective  and  arbitrary  nature  of  these  notes  makes  effective  follow-up  exams 
problematic. 

We  have  conducted  a  study  to  advance  fundamental  understanding  of  palpation  and  solve  these  practical 
problems  through  the  creation  of  a  new  tactile  mapping  device  (TMD).  This  device  will  measure  three  key 
variables  during  palpation:  the  examiner’s  search  patterns,  the  applied  forces,  and  the  small-scale  pressure 
variations  at  the  skin  due  to  lumps.  We  have  integrated  and  pilot-tested  a  prototype  TMD  resisting  0f  a 
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novel  three-dimensional  (3-D)  camera  that  can  track  finger  motion  and  breast  deformation  in  video  speed, 
a  six  degree-of-freedom  (DOF)  force/torque  sensing  device  measuring  the  applied  forces,  and  a  novel  pressure 
distribution  sensor  array  that  can  image  tactile  profile  of  lumps.1  The  primary  objective  aims  that  (1)  new  tactile 
mapping  technology  can  quantitatively  measure  the  location  and  applied  forces  in  breast  palpation,  and  the 
tactile  features  of  detected  breast  lumps;  and  (2)  new  device  can  accurately  characterize  and  document  breast 
lumps  and  will  improve  clinicians’  ability  to  monitor  changes  in  lump  across  time  and  possibly  to  distinguish 
mahgnant  from  benign  lumps.  From  a  set  of  “images”  of  the  suspect  mass,  a  neural  network  supported  pattern 
analysis  system  will  extract  the  invariant  properties  of  the  lump,  such  as  the  depth  and  size,  based  on  a  nonlinear 
model  of  sensor-tissue  interaction  with  hard  inclusions.  This  TMD  system  will  make  it  possible  for  the  first 
time  to  quantitatively  and  objectively  record  and  characterize  the  processes  and  findings  of  breast  palpation. 
While  initially  the  system  will  be  used  to  perform  clinical  breast-  examination,  we  believe  that  eventually  it  may 
be  used  in  breast  self-examination  by  women  through  tele-home  care.4,5 

2.  SYSTEM  OVERVIEW 

The  goal  of  this  research  is  to  extend  the  range  and  resolution  of  breast  palpation  methods,  thus  increasing 
palpation  sensitivity  and  specificity  (i.e,  the  ability  to  detect  lumps  and  distinguish  clinically  significant  changes). 
Using  3-D  camera^  force  sensor,  and  tactile  arrays,  we  can  create  reproducible  tactile  maps  (or  images)  of  the 
palpable  abnormalities  of  the  breast.  The  novel  3-D  camera  is  able  to  provide  real-time  3-D  motion  measurement 
(for  tactile  probe)  and  surface  profile  measurement  (for  breast  tissue  deformation)  at  a  rate  of  30  frames  per 
second.  We  have  further  integrated  and  evaluated  related  sensing  technologies  including  tactile  sensor  array 
and  six  DOF  force/torque  sensor  to  extract  tactile  information  from  breast  palpation,  where  the  sensor-tissue 
interaction  of  detected  palpable  lesion  can  be  quantitatively  measured  and  displayed  in  terms  of  tactile  images. 
Our  laboratory  experiments  aimed  at  characterizing  the  simulated  tumors  in  breast  phantom  have  demonstrated 
that  palpable  inclusions  can  be  located  to  within  1  mm 

Our  prototype  system  incorporates  a  3-D  Rainbow  camera,  a  tactile  sensor  array  probe,  a  breast  model  with 
lranPs>  a  force/torque  sensor  on  which  the  breast  model  is  mounted,  and  a  graphical  user  interface 
(GUI).  Figure  1  shows  a  block  diagram  of  the  tactile  mapping  device  system  to  be  used  for  measurement  of 
palpable  abnormalities  in  the  laboratory.  In  a  typical  task,  upon  detecting  a  suspicious  lump,  the  examiner  will 
bring  the  TSAP  into  contact  with  the  tissue  at  the  palpation  site.  For  a  thorough  tactile  mapping  procedure  this 
involves  3-D  positioning  the  probe  through  the  camera  facing  the  site.  The  resulting  pressure  distribution  a^oss 
the  probe  surface  is  measured  by  a  Tactile  Sensor  Array  Probe  (TSAP)  with  associated  readout  electronics. 
Multiple  tactile  images  with  various  force  levels  and  torque  angles  will  be  required  for  the  lesion  characterization. 

A  computer  processes  the  signals  to  generate  appropriate  output  for  visual  display  on  a  monitor  or  raw  data  to 
the  physician’s  office  through  a  telecommunication  channel.  Below,  we  describe  the  components  of  the  tactile 
mapping  system. 

2.1.  3-D  Rainbow  Camera 

For  the  accurate  positioning  of  the  TSAP,  a  novel  3-D  rainbow  camera  is  adopted  in  our  system,  which  is  suitable 
for  this  high-speed  3-D  machine  vision  application.  Figure  2  shows  the  rainbow  camera  and  the  acquired  sample 
range  image.  It  exploits  a  color  light  projector  to  illuminate  the  object  in  the  scene  and  using  an  off-the-shelve 
color  camera  to  obtain  a  full-frame  color  image  of  the  scene.  The  color  of  the  projecting  light  with  spatially 
continuously  varying  wavelength  is  encoded  with  information  of  the  corresponding  projection  angle.  Each  pixel 
of  the  color  image  is  associated  with  a  unique  ray  through  the  focal  point  of  the  camera.  Since  the  angle 
between  camera  axis  and  the  ray  is  known  the  resulting  angle-side-angle  triangulation  problem  can  be  easily 
solved- when  the  distance  between  the  light  projector  and  the  camera  is  fixed;  Thus,  using  only  one  camera, 
the  full  frames  of  3-D  range  images  can  be  obtained  directly  at  the  camera  frame  rate.  The  spatially  varying 
wavelength  light  is  generated  by  a  white  light  passing  through  a  linear  variable  wavelength  filter.  We  have 
performed  experiments  to  investigate  the  actual  range  accuracy  and  the  major  error  contributors.  The  results 
show  that  the  3-D  profile  of  test  artifacts  were  less  than  1  mm.  The  spatial  resolution  of  the  system  is  limited 
only  by  the  spatial  resolution  of  the  camera  optical  sensing  element  and  is  currently  able  to  provide  at  least 
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Figure  3:  A  drawing  of  a  tactile  array  sensing  (adapted  from  Son  1996). 

1024x1024  pixel  resolution.  When  a  palpation  site  is  determined,  several  3-D  range  images  will  be  acquired  to 
record  the  site  information. 

2.2.  Tactile  Sensor  Array  Probe  (TSAP) 

The  core  component  of  the  TMD  is  a  highly  sensitive  tactile  sensor  array  probe.  Figure  3  illustrates  the 
operation  of  the  capacitive  tactile  array  sensor  used  in  this  system.  It  shows  a  drawing  of  the  array  with  the 
copper  layers  and  the  silicone  rubber  spacers,  which  is  based  on  the  design  by  Fearing6  and  Howe.4  The  array 
is  composed  of  two  crossed  layers  of  copper  strips  separated  by  thin  strips  of  silicone  rubber.6  Each  crossing  = 
area  forms  a  capacitor,  and  when  a  force  is  applied  onto  where  the  strips  cross,  the  distance  between'  the  strips 
decreases  and  the  capacitance  increases.7  Specially  designed  electronics  will  measure  the  capacitance  of  each 
element  and  relate  the  capacitance  change  to  the  force  applied  to  each  element.  By  measuring  the  capacitance 
variations  from  all  the  elements  simultaneously,  we  can  determine  the  spatial  distribution  of  pressure  across  the  ' 
sensor  array.  .  -4*  •  ’  '  V  ' 

Figure  4  shows  a  photograph  of  a  commercial  prototype  TSAP.  The  sensor  Array  is  made  with  ah  inexpensive, 
photolithography  and  etching  process  and  can  be  easily  attached  to  a  variety  of  probe  shapes.  In  this  prototype 
specification,  it  is  composed  of  an  eight  by  eight  tactile  sensor  with  elements  that  are  4  mm  on  a  side.  The 
sensor  is  mounted  on  a  plastic  brass  backing  plate  with  a  surface  that  has  been  machined  into  a  section  of  a 
square.  The  backing  plate  is  5.08  cm  on  a  side  and  the  effective  sensing  area  is  3.20  cm  on  a  side.  We  decided 
to  make  the  sensor  flat  in  order  to  minimize  inhomogeneity  because  the  resulting  pressure  distribution  should 
have  a  uniform  overall  signal  to  noise  ratio.  The  tactile  images  will  then  be  consistent  when  used  for  various 
breast/chest  background  textures.  The  spatial  resolution  of  the  tactile  array  is  4  mm  long  where  the  smallest 
masses  that  we  are  currently  interested  to  characterize  are  on  the  order  of  1  cm  in  diameter.  Smaller  elements 
would  increase  spatial  resolution  at  the  cost  of  lower  coverage  area  and  low  sensitivity  since  the  capacitance  is 
proportional  to  the  element  area.7  The  tactile  sensor  has  been  shielded  to  provide  electronic  insulation- 

2.3.  Force/Torque  Sensing  System 

The  relationship  between  the  hard  inclusion  (i.e.,  lump)  and  the  perdeved  tactile  image  from  the  TSAP  is 
nonlinear  and  complex.  In  order  to  characterize  and  later  extract  the  tactile  features  of  the  detected  breast 
lumps,  the  TMD  operation  requires  that  the  forces/torques  exerted  by  the  operator  of  the  TSAP  on  the  breast 
model  be  measured  together  with  the  corresponding  tactile  images.  It  would  be  ideal  if  a  force/torque  sensor 
were  to  be  mounted  between  the  wrist  and  the  hand  of  the  operator.  Since  it  is  problematic  and  impractical 
to  mount  a  force/torque  sensor  in  such  a  configuration,  we  dedded  to  mount  it  under  a  base  which  the  breast 
model  is  placed  on.  The  forces/torques  acquired  by  the  sensor  in  this  configuration  can  be  transformed  to  those 
exerted  by  the  operator  via  proper  coordinate  transformation.  Figure  5  shows  a  breast  model  laying  on  a  round 
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Figure  5:  JRS  Force/Torque  Sensor. 


tese  under  winch  a  JR?  force/torque  is  mounted.  The  JR3  force/torque  mainly  consists  of  a  JR3  monolithic 
six  DOF  force  sensor  and  a  JR3  Intelligent  Support  System,  Comprised  of  boards  for  signal  conditioning,  data 
acquisition  and  processing.  A  computer  program  was  written  in  Visual  C++  to  request  the  JR3  force/torque 
sensing  system  to  send  to  the  computer  six  forces/torques  along  about  the  Cartesian  x-,  j-,  z-zx.es  assigned  to 
the  base.  The  transmitting  and  receiving  of  data  between  the  computer  and  the  force/torque  system  is  carried 
out  through  serial  port. 

2.4.  Breast  Models 

The  breast  models  used  in  our  laboratory  experiments  were  provided  by  the  HEALTH  EDCO,  a  Division  of 
WRS  Group,  Inc.  The  models  were  made  from  BIOLIKE™  synthetic  tissue  that  feels  just  like  a  real  breast 
[IJ.  Two  models  were  used  in  our  tests  to  collect  preliminary  data.  Each  of  them  has  5  lumps  that  simulate 
easy-  and  hard-to-find  breast  tumors  with  various  sizes  and  depths  at  different  locations.  One  of  the  models 
is  the  geriatric  breast  that  is  ideally  shaped  to  address  the  special  problems  of  older  women.  The  geriatric 
model  simulates  the  natural  stretching  of  the  tissue  with  age.  Figure  6  show  the  distribution  of  simulated  lumps 
(lesions)  in  the  breast  model  and  the  historical  spatial  distribution  of  the  breast  cancer  respectively;  It  should 
be  noticed  that  as  high  as  50%  of  breast  cancers  occur  at  the  axilla  area  which  is  difficult  to  imaging, 

2.5.  Graphical  User  Interface 

For  both  tactile  documentation  and  interactive  training,  a  graphical  user  interface  (GUI)  plays  a  very-important 
role  of  optimizing  both  components  and  examiner’s  performance.  In  our  prototype,  Microsoft  Visual  0++  was 
employed  to  implement  the  GUI  to  visually  display  pertinent  data  of  the  TSAP  and  the  force/torque  sensor. 
Figure  7  shows  the  output  window  of  the  GUI  that  presents  on-line  the  pressure  distribution  acquired  by  the 
TSAP,  the  time  response  of  the  6-DOF  force/torque  apphed  by  the  examiner.  The  user  can  interact  with  the- 
TMD  functions  through  the  designed  menu 


of  the 


of  breast  cancer 


ce/torque  feedback  (upper  left),  tactile 
:  left),  and  3D  mesh  contour  map  based 
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3.  EXPERIMENTAL  RESULTS 

!°  danof'ate  ““  eSfcto»«a  of  a»  TMD  system  for  the  characterisation  and  documentation  of 
detected  lumps,  intensive  experiments  have  been  conducted  to  test  tbe  sensitivity  and  reproducibility  of  the 
SSL"  meamirmg  the  applied  farces  and  pressure  profiles  due  to  lumps.  The  priuuuy  ob^Stee 1st 
^erform^i^  rf“h^^e  of  P«^  Wans  using  the  TMD  system  and  to  evaluate  the 

Dalrato^i^'s™”™',5  Conduf!d  ,Mot’  ;hc  forces/torques  applied  by  the  human  examiner  during  the 
palpation  are  recorded  m  terms  of  six  force  parameter  -rallies,  i.e  Fx  Fv  Fz  Mx  Mv  M?  Pnr  r.8,., 

imaging  of  simulate!  lesions,  the  TSAP  is  used  to  acquire  various  tkSr  £££ 

h? Z**  *  ““  W  ^  «»  rdationshi^  bSSftSS. 

charactenstics  are  expected  to  be  complex  and  nonlinear,  we  believe  that  the  inverse  problem  (extraction 

r fe^Wrtf3 7  ^  t^tile  ™ages)  0811  be  8oived  when  ^cient  tactile  information  is  acquired  As  a 
result,  for  each  of  the  lesions  after  pushing  the  TSAP  against  the  breast  to  achieve  certain  level  offeree  the  ex 

amrner  rotates  the  TSAP  to  five  different  orientations  at  each  of  which  the  information  about  the  forces/toroues 
and  tactile  image  axe  simultaneously  recorded.  OKI  is  the  initial  orientation  of  the  TSAP  after  being  pushed  by 

fo*TT7ZZ!°  TiT  le7l0fJ0ICe-  °m>  0R3’  0R4  811(1  0R5  ««  TSAP  orientations  after  it  k  routed 
forward,  backward*  left  and  right,  respectively.  ^ 

n»me ^aXT7  oUlthe  f?°!lprocedure  on  t"°  lesions:  Lesion  1  and  Lesion  4  as  identified  in  Ffeure  6 

meJur^dbvtheTSAP?nT  (!?ft)  “d the  shadin°  ^erpolation  (right)  of  the  pressure  distribSion 

q-dz  Q  ,  ,  r  tlie  teP  r0w  for  OR1  and  the  following  rows  for  OR2,  OR 3,  OR4  and 

OR5.  Figure  9  shews  the  results  obtsiued  for  Lead  4  presented  in  the  same  manner  as  ItaiCi 
of  farces/torques  for  the  experiments  conducted  for  Lesions  1  and  4  are  tabulated  in  Table 

Table  1:  Results  force/ torque  measurement  for  Lesion  1  and  4, 


Lesion  TSAP  Orientation 


6DOF  Force  fe  Torque  (Ib./imlb.) 


OR1 

Fa 

Fy 

F  . 

M* 

Mj, 

M* 

-0.89 

0.44 

•43.13 

-2.58 

-10.80 

0.50 

OR2 

-2.34 

2.50 

-0.10 

-5.26 

-22.02 

1.59 

OR3 

1.69 

2.96 

-0.03 

-2.28 

-10.50 

3.57 

OR4 

3.20 

-0.38 

0.27 

0.40 

0.40 

2.09 

OR5 

-1.83 

0.38 

-0.03 

-2.38 

-9.32 

1.99 

OR! 

-3.17 

1.77 

-0.16 

5.58 

25.59 

-0.89 

OR2 

-2.08 

2.48 

-0.22 

4.37 

18.35 

1.78 

OR3 

-8.81 

-0.03 

-0.09 

1.98 

10.71 

0.20 

OR4 

-3.67 

-1.72 

0.00 

2.18 

.10.01 

0.39 

OR5 

-4.77 

4.02 

-0.34 

3.18 

10.81 

3.97 

at  tte  expemnent  ”  ^0nJ1  ^  4-  «  apply  «*  same  procedure  ou  a  location 

Firnreg  n0  lesron  The  resulting  display  shows  no  peat.  Examining  the  graphical  displays  in 

Fjure  8  and  9  in .comparison  with  that  of  the  case  of  no  lesion,  we  notice  that  a  peat  comSiaJo^STth” 

S  Won  Pf™y  “f-ion  iere  that  a  peat  in  the  graphical  display  indicates  SSsZL 

f  force  and  torque  that  we  recorded  in  conjunction  with  the  tactile  sensing  during 
°n  “  46  l0Cati0n  “d  iC”  “a*  we  apply  on  L  harl 


4.  DISCUSSION  AND  CONCLUSION 

““  deteCti0‘1  °fiumps  fa  ^ast  Palpation  using  a  tactile  mapping  device  (TMD) 
system.  It  presented  the  mam  components  of  the  TMD  and  explained  their  operations.  B^eriments  ware 
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Figure  8:  Display  of  experimental  results  for  lesion  1. 


conducted  on  a  breast  model  to  detect  two  different  types  of  lumps.  Experimental  data  showed  that  the  peat  in 
the  display  of  pressure  distribution  of  the  TMD  indicated  the  presence  of  a  lump.  The  results  demonstrated  the 
feasibility  of  the  TMD  for  improving  breast  examination  technique  in  diagnosis,  documentation,  and  training.  In 
particular,  the  results  have  shown  that  new  tactile  mapping  technology  can  quantitatively  measure  the  location 
and  applied  forces  in  breast  lumps;  the  prototype  interactive  training  program  can  track  finger  motions  and 
applied  forces  during  breast  palpation  in  which  on-line  feedback  can  help  the  training  to  better  understand  the 
search  strategy  and  adjust  force  level  to  increase  the  sensitivity.  "With  the  proven  power  of  nonlinear  signal 
processing  using  both  convolution  neural  networks  (CNN)  and  multi-layer  perception  (MLP),  we  expect  that 
the  tactile  parameters  associated  with  the  lesions  (i-e.,  the  size  and  depth  of  the  lesion)  csn  be  estimated  more 
accurately  than  those  by  the  conventional  approaches  using  the  first  principle  of  engineering.7 
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ABSTRACT 

Non-rigid  image  registration  is  a  prerequisite  for  many  medical  imaging  applications  such  as  change  analysis  in 
image-based  diagnosis  and  therapy  assessment.  Nonlinear  interpolation  methods  may  be  used  to  recover  the 
deformation  if  the  correspondence  of  the  extracted  feature  points  is  available.  However,  it  may  be  very  difficult 
to  establish  such  correspondence  at  an  initial  stage  when  confronted  with  large  and  complex  deformation,  In 
this  paper,  a  mixture  of  principal  axes  registration  (mPAR)  is  proposed  to  tackle  the  correspondence  problem 
using  a  neural  computation  method.  The  feature  is  to  align  two  point  sets  without  needing  to  establish  the 
explicit  point  correspondence.-  The  mPAR  aligns  two  point  sets  by  minimizing  the  relative  entropy  between 
their  probability  distributions  resulting  in  a  maximum  likelihood  estimate  of  the  transformation  mixture.  The 
neural  computation  for  the  mPAR  is  developed  using  a  committee  machine  to  obtain  a  mixture  of  piece-wise 
rigid  registrations.  The  complete  registration  process  consists  qf  two  steps:  (1)  using  the  mPAR,  to  establish 
an  improved  point  correspondence  and  (2)  using  a  multilayer  perceptron  (MLP)  neural  network  to  recover 
the  nonlinear  deformation.  The  mPAR  method  has  been  applied  to  register  a  contrast-enhanced  magnetic 
resonance  (MR)  image  sequece.  The  experimental  results  show  that  our  method  not  only  improves  the  point 
correspondence  but  also  results  in  a  desirable  error-resilience  property  for  control  point  selection  errors. 

Keywords:  M  ixture  of  Principal  Axes  Registration,  Non-Rigid  Image  Registration,  Change  Analysis,  Neural 
Computation,  Multilayer  Perceptron,  Control  Points  Selection. 

1.  INTRODUCTION 

Image  registration  is  an  essential  step  for  many  medical  image  analysis  applications  such  as  image  fusion,  quan¬ 
titative  change  analysis  in  image-based  diagnosis  and  therapy  assessment.1’2  For  example,  medical  diagn^c 
often  benefits  from  the  complementary  information  in  images  of  different  modalities.  Multimodality  ima&ng  of 
breast  is  regarded  as  a  powerful  diagnostic  tool.  In  breast  imaging,  different  modalities  offer  different  diagnostic 
information.  In  order  to  fuse  the  complementary  information,  the  alignment  of  the  images  is  the  first  step  to 
be  performed  to  overcome  the  complication  of  misregistrations  caused  by  patient  motion  and  physical  change. 
Image  registration  consists  of  aligning  the  images  by  scaling,  rotating,  and  translating,  one  or  both  images  so 
that  they  axe  of  the  same  size  and  have  the  same  orientation  and  location.  Mathematically,  image  registration 
can  be  formulated  as  to  estimate  the  transformational  geometry  from  two  feature  point  sets,  i.e.,  to  recover  a 
matrix  representation  requiring  a  set  of  correspondence  matches  between  features  in  the  two  coordinate  sys¬ 
tem.  Arun  et.  al.  present  an  algorithm  for  finding  the  least-squares  solution  of  the  transformation  matrix, 
which  is  based  on  the  decoupling  of  translation  and  rotation  and  the  singular  value  decomposition  (SVD)  of  a 
cross-covariance  matrix. 
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The  major  limitation  of  the  above  method  is  twofold:  (1)  while  feature  matching  methods  can  give  quite 
accurate  solutions,  obtaining  correct  correspondences  of  features  is  a  hard  problem,  especially  in  the  case  of 
images  acquired  using  different  modalities  or  from  inter-subjects;  and  (2)  a  rigidity  assumption  is  heuristically 
imposed  that  is  usually  not  valid  in  medical  images  with  non-rigid  human  organs  such  as  breast.2’4  One 
popular  method  that  does  not  require  feature  correspondence  is  the  principal  axes  registration  (PAR)  method.1 
The  PAR  method  is  based  on  the  relatively  stable  geometric  properties  of  image  features,  i.e.,  the  geometric 
information  contained  in  these  stable  image  features  is  often  sufficient  to  determine  the  transformation  between 
images.2  However,  the  PAR  method  results  in  a  rigid  body  transform  with  scaling  coefficients.  The  transform 
works  effectively  if  the  objects  to  be  registered  have  orientation  differences  within  a  certain  limit,  dependent  on 
the  objects >  shape  and  if  their  shapes  axe  well  matched.  If  their  orientations  axe  too  different,  the  transform 
may  arrive  at  the  wrong  rotation  angle  due  to  problems  with  object  symmetry.  Once  again,  the  method  cannot 
handle  the  cases  with  non-rigid  objects. 

in  this  paper,  we  present  a  neural  computation  based  non-rigid  image  registration  method  using  piece¬ 
wise  rigid  transformation.  In  other  words,  rather  than  using  a  single  transformation  matrix  that  gives  rise 
to  a  large  registration  error,  we  attempt  to  interpolatively  apply  a  mixture  of  transformations  to  minimize  the 
registration  error.  By  generalizing  PAR  to  a  mixture  of  principal  axes  registration  (mPAR)  scheme,  the  mixture 
is  fit  using  the  expectation-maximization  (EM)  algorithm  by  performing  a  soft  partitioning  of  the  data  set.  The 
registration  process  consists  of  two  steps:  (1)  using  the  mPAR  to  establish  an  improved  point  correspondence 
and  (2)  using  a  multilayer  perception  (MLP)  neural  network  to  recover  nonlinear  deformation  based  on  the 
established  point  correspondence.  A  finite  mixture  of  transformations  is  obtained  through  aligning  the  two 
point  sets  by  minimizing  the  relative  entropy  between  their  probability  distributions  resulting  in  a  mayirrmm 
likelihood  estimate  of  the  transformation  matrix.  Specifically,  a  probabilistic  adaptive  principal  components 
extraction  (PAPEX)  algorithm10  is  developed  to  estimate  the  transformation  matrix  using  the  orthogonal  set 
of  eigenvalues  and  eigenvectors  of  the  auto  covariance  matrix.  A  committee  machine  is  then  used  to  combine 
multiple  transformations  for  the  recovery  of  the  total  transformational  geometry  of  the  non-rigid  object.  An 
accurate  point  correspondence  can  be  established  by  applying  the  finite  mixture  of  transformations  to  the  two 
point  sets.  Finally,  an  MLP  neural  network  is  adopted  to  recover  the  nonlinear  deformation  by  the  polynomial 
mapping  function  based  on  the  point  correspondence.  '  • 

The  mPAR  method  has  been  applied  to  register  a  contrast-enhanced  breast  magnetic  resonance  image  (MRI) 
sequence.  With  MRI,  the  breast  is  usually  imaged  functionally  using  a  contrast  agent.  The  rate  of  uptake  of 
the  contrast  agent  is  used  to  characterize  a  given  breast  tissue.  This  rate  is  studied  by  observing  the  change  in 
gray  level  of  the  tissue  at  several  instants  of  time  after  injection  of  the  agent.  However,  the  non-rigid  nature 
of  breast  tissue  almost  guarantees  that  these  different  images  would  be  misregistered.  To  study  the  uptake 
process  effectively  and  to  also  make  a  quantitative  estimate  of  change  in  uptake,  it  would  be  necessary  to  align 
or  register  these  different  images.  Image  registration  is  thus  an  important  problem.  To  align  these  different 
images,  we  first  extract  the  control  objects  of  the  image  using  a  stochastic  segmentation  method  With  the 
extracted  control  objects,  we  then  apply  our  mPAR  method  to  recover  the  total  transformational  geometry 
of  the  contrast-enhanced  MRI.  An  accurate  point  correspondence  can  be  established  by  applying  the  mixture 
of  transformations  to  the  two  point  sets.  As  a  final  step,  we  apply  an  MLP  neural  network  to  recover  the 
nonlinear  deformation  in  the  form  of  the  polynomial  transformation  using  the  established  point  correspondence. 
The  experimental  results  shows  that  the  image  differencing  after  registering  two  images  by  applying  our  mPAR 
method  can  greatly  improve  the  accuracy  in  extraction  of  the  enhanced  area  compared  to  that  of  without  any 
registration. 


2.  THEORY  AND  METHOD 

Assume  two  3-D  data  point  sets  {pij4}  and  {piS},  i  =  1, 2, AT,  are  related  by 


PiB  =  RPiA  +  T  +  Ni  (1) 

where  R  is  a  rotation  matrix,  T  a  translation  vector,  and  N *  a  noise  vector.  Given  {piA}  and  {piB},  Arun 
et  d.  present  an  algorithm  for  finding  the  least-squares  solution  of  R  and  T  by  decoupling  translation  and 
.rotation  and  using  the  singular  value  decomposition  (SVD)  of  a  3  x  3  cross-covariance  matrix.3 
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Suggested  by  information  theory,5  we  can  consider  the  feature  point  sets  in  two  images  as  two  separate 
realizations  of  the  same  random  source.  Hence,  to  align  two  sets  of  feature  points,  we  do  not  need  to  establish 
point  correspondences  to  extract  the  transformation  matrix.  Instead,  we  can  aligning  the  two  point  sets  by 
minimizing  the  relative  entropy  between  their  probability  distributions.  In  other  words,  if  P{Pi}  denotes  the 
distribution  of  the  feature  point  set  in  an  image,  we  have  the  simple  probabilistic  relationship  to  describe  the 
transformation  between  two  point  sets: 

F{Pis)  =  P{RPiJ,+T}  +  e  (2) 

where  e  is  the  noise  component.  Since  the  probability  distributions  can  be  computed  independently  on  each 
image  without  the  need  to  establish  feature  correspondences,  and  given  the  two  distributions  of  the  control 
point  sets  in  the  two  images,  we  can  recover  the  transformation  matrix  in  a  simple  fashion,2  as  we  now  sketch. 

2.1.  Registration 

For  observation  of  the  distributions,  we  can  estimate  R  and  T  by  minimising  the  relative  entropy  (Kullback- 
Leibler  distance)  between  P{PjBy  and  P{ rp<a+t}-  The  least  relative  entropy  estimator  is  defined  as 

arg^D(P{p.B}|[P{RpiA+T})  (3) 

where  P  denotes  the  relative  entropy  measure.  Following  the  same  strategy  to  decouple  translation  and  rotation 
as  in,  _  we  can  define  a  new  data  point  by  qlA  =  piA  -  P°  and  q,B  =  PjB  -  p|,  where  p^  and  p|  are  the 
centroids  of  {p^}  and  {P;b},  respectively.  Then  the  optimal  geometric  transformations,  R  and  T,  can  be 
computed  as 


R  =  USHU^ 
T  =  p|-Rp° 


(4) 

(5) 


where  the  superscript  t  denotes  matrix  transposition,  XJA  and  U %  axe  3x3  orthonormal  matrices,  and  H 
is  a  3  x  3  diagonal  matrix  with  element  hm  ==  \/A Note  that  the  transformation  U  consists  of  the 
orthonormal  set  of  eigenvectors  and  hm  is  the  squared  root  of  the  eigenvalues  Am  of  the  auto-covariance  matrix 
C  for  m  =  1, 2, 3  and  for  {p^}  and  respectively.  ‘ 

However,  because  of  its  global  linearity,  the  application  of  PAR  is  somewhat  limited  to  deal  with  rigid 
objects  only.  An  alternative  paradigm  is  to  model  a  multimodal  control  point  set  with  a  collection  of  local 
linear  models.  The  method  is  a  two-stage  procedure:  a  soft  partitioning  of  the  data  set  followed  by  estimation 
of  the  principal  axes  within  each  partition.7  Recently  there  has  been  considerable  success  in  using  standard 
finite  normal  mixture  (SFNM)  to  model  the  distribution  of  a  multimodal  data  set.  The  association  of  a  SFNM 
distribution  with  PAR  offers  the  possibility  of  being  able  to  register  two  images  through  a  mixture  of  probabilistic 
principal  axes  transformations. 7  Assume  that  there  are  Kq  control  point  clusters,  where  each  control  point 
cluster  defines  a  transformation  {R/p,T&}.  For  a  point  pnA,  its  new  locations  corresponding  to  each  of  the 

transformations  are  Pn&  =  RfcPnA  +  T&  for  k  =  1, Kq.  Further  assume  that  the  control  point  set  defines  a 
SFNM  distribution 

z* 

/(pi) = cfc)  (6) 

s.  k=1 

where  g  is  the  Gaussian  kernel  with  mean  vector  (ik  and  auto-covariance  matrix  Cfc,  and  Wfc  is  the  mixing  factor 
proportional  to  the  number  of  control  points  in  cluster  k.  For  each  of  the  control  point  sets  {pm}  and  {pig}, 
the  mixture  distribution  is  fit  using  the  expectation-maximization  (EM)  algorithm.  In  principle,  the  E  step 
involves  assigning  to  the  linear  models  responsibilities,  i.e.,  the  posterior  Bayesian  probability,  from  the  control 
points;  the  M  step  involves  re-estimating  the  parameters  of  the  linear  models  from  the  above  assignment.7 

Thus  the  statistical  membership  of  point  Pn^  belonging  to  each  of  the  control  point  dusters  can  be  derived 
by 

Znfc  =  P(Rfc.Tfclp^)  =  l^iPnA^C^)  (7) 
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the  posterior  probability  of  {R*,  T*}  given  p nA-  We  can  define  the  mPAR  transformation  as 


K0 

Pn  =  ZnkPnk 

k=l 

=  v  c«)  (R,p^ + Tt) 

fc=l  f\PnA) 


(9) 


where  {R*,  T*J  is  determined  based  on  Cfc>i}and  {^*3,  C^b},  which  we  have  estimated  in  the  previous 
step  using  the  EM  algorithm.  This  philosophy  for  recovering  transformational  geometry  of  the  non-rigid  objects 
is  similar  in  spirit  to  the  divide- and- conquer  principle,6  under  which  the  relative  entropy  between  the  two  point 
sets  reaches  its  minimum 

arg  D(P{pjB}\\P{E^  ^(RtPM+Tt)})  (10) 

both  globally  and  locally. 

Based  on  a  mixture  of  probabilistic  principal  axes  transformations,  the  next  section  describes  a  neural 
computation  using  a  committee  machine  to  obtain  a  mixture  of  piece-wise  rigid  registrations,  which  gives  a 
reliable  point  correspondence  using  multiple  extracted  object  control  points. 


2.2.  Neural  Computation 

There  are  many  numerical  techniques  to  perform  the  maximum  likelihood  estimation.  The  most  popular  method 
is  the  EM  algorithm.  The  EM  algorithm  first  calculates  the  posterior  Bayesian  probabilities  of  the  data  through 
the  data  observations  and  the  current  parameter  estimates  (E-Step)  and  then  updates  parameter  estimates 
(M-Step).  The  procedure  cycles  back  and  forth  between  these  two  steps.  A  neural  network  interpretation  of 
the  EM  algorithm  is  given  in.8  Because  of  its  reputation  of  being  slow  in  which  new  information  acquired  in 
the  expectation  step  is  not  used  immediately,  on-line  versions  of  the  EM  algorithm  are  proposed  for  large-scale 
sequential  learning.  We  adopt  a  fully  unsupervised  and  incremental  stochastic  learning  algorithm  to  implement 
the  EM  algorithm.  The  scheme'  provides  winner-takes-in  probability  (Bayesian  “soft”)  splits  of  the  control 
points,  hence  allowing  the  data  to  contribute  simultaneously  to  multiple  clusters.  The  incremental  stochastic 
learning  EM  algorithm  can  be  described  as  follows: 


E-Step 

»l')g(pi-nl<4‘>,ci.‘)) 

M  —  Step 

■  Vk+1)  =  A 4°  +  o(0(p»+l  - Mfcl))z(i+l)fc> 

C<<+1)  =  c«  +  6(i)[(pi+1  -  /4°)(Pi+i  -  /4°)T  -  cfhi+D*/ 


(11) 


(12) 

(13) 


Tk 


<i+1>  - 


i  +  1 


+  7TT^i+1^ 


(14) 


for  k  =  1, Kq  and  for  {p^}  and  {pis},  respectively.  In  the  above  equations,  a(i)  and  6(i)  are  introduced 
as  the  learning  rates  that  converge  to  zero,  ensuring  unbiased  estimates  after  convergence.  This  procedure  is 
termed  as  neural  computation  of  the  EM  algorithm.  At  each  complete  cycle  of  the  algorithm,  we  first  use 
“old”  set  of  parameter  values  to  determine  the  posterior  probabilities  Z(i+i)fc*  These  posterior  probabilities  are 

then  used  to  obtain  “new”  values  7t£+1\  The  algorithm  cycles  back  and  forth  until  the  value  of 

relative  entropy  between  the  data  histogram  and  mixture  model  reaches  its  minimum 


arg  min  D(P{p.}||/(pi)) 


(15) 
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for  {pia}  and  {pis}s  respectively. 

With  a  soft  partitioning  of  the  data  set  using  Eqs.  (11-14),  control  points  mil  now  effectively  belong  to 
more  than  one  cluster  spatially.  Thus,  the  effective  input  values  are  p &  —  Z&  (p4  -  pk)  for  an  independent 
registration  transformation  k  in  the  committee  machine.9  We  then  extend  an  adaptive  principal  components 
extraction  (APEX)  algorithm  to  a  probabilistic  version,  Le.,  PAPEX,10  to  determine  Ufc  for  {piA}  and  {piB}, 
respectively,  summarized  as  follows.  *  5 

1.  Initialize  the  feedforward  weight  vector  u*,*  for  m  =  1, 2, 3,  and  the  feedback  weight  vector  a**  to  small 
random  values  for  m  2, 3,  at  time  i  =  1.  Assign  a  small  positive  value  to  the  learning  rate  parameter  i j. 


2.  Set  m  =  1,  and  for  %  =  1,2, ...,  compute 

yu(i)  =  «u(Ozife(Pi~Mfc)  (16) 

uu(t  + 1)  =  ulfc(i)  +  ■n[ylk(i)zlk(pi  -  fik)  -  ylk(i)nlk(i)}  (17) 

For  large  i  we  have  Ujfc(i)  >  where  uu,  is  the  eigenvector  associated  with  the  largest  eigenvalue  of 
the  cluster  k,  and  Xlk  =  jj  tfk(i). 

3.  Set  m  =  2,  and  fox  %  =  lf  2}  compute 

y(m-l)k(i)  =  [m(i), &*(*),.. .,y(m-i)k(i)]T  (18) 

Wmfc(*)  =  U^k(i)zik(pi  —  fik)  -I-  a^j,(i)y^m_1)k(i)  (19) 

Um*(i  +  l)  =  umS;(i)-|-)7[ymfc(j)2ifc(pi-/LtA.)-t4Jfc(i)uTnjfc(i))  (20) 

a mk(i  +  1)  =  -  Tl[ymk(i)y(m-i)k(i)  +  2/^(i)amA(i)]  (21) 


For  large  i  we  have  U2k(i)  >  1x2*,  where  112*  is  the  eigenvector  associated  with  the  second  largest 
eigenvalue  of  the  cluster  k,  and  X2k  =  j?  vUS)- 

4.  Set  m  =  3,  go  to  step  3.  For  large  i  we  have  o.zk(i)  — *• 113^.  where  1134  is  the  eigenvector  associated  with 
the  third  largest  eigenvalue  of  the  cluster  k,  and  X$k  =  2i=i  V% *(*)•  The  next  step  is  to  introduce  a 

committee  machine  for  combining  multiple  transformations  in  a  mixture  fashion. 

Consider  the  modular  network  nature  of  Eq.  (8),  we  can  easily  realize  that  it  is  a  mixture  of  local  expert 
model.  Each  local  expert  performs  supervised  learning  through  our  PAPEX  algorithm,  and  the  individual 
responses  of  the  experts  are  nonlinearly  combined  by  means  of  a  single  gating  network.  This  is  our  committee 
machine  to  combine  multiple  transformations  to  form  the  total  transformational  geometry  of  the  non-rigid 
object.  The  neural  computation  of  a  committee  machine  can  be  achieved  by  distributing  the  teaming  tasks, 
among  a  number  of  experts,  which  in  turn  partitioning  the  input  space  into  a  set  of  subspaces.  The  experts 
are  in  theory  performing  supervised  learning  in  that  the  individual  outputs  are  combined  to  model  the  desired 
response.  There  is,  however,  a  sense  in  which  the  experts  are  also  performing  self-organized  learning;  that  is, 
they  self  organize  to  find  a  good  partitioning  of  the  input  space  so  that  each  expert  does  well  at  modeling  its  own 
subspace,  and  as  a  whole  group  they  model  the  input  space  jointly.  It  is  assume  here  that  the  different  experts 
work  best  in  different  regions  of  the  input  space  in  accordance  with  the  probabilistic  generative  model.  In  our 
case,  the  effective  input  values  are  p-k  =  £&(p*  —  fik)  fox  an  independent  registration  transformation  k  in  the 
networks.  A  committee  machine  consists  of  k  supervised  modules,  with  a  soft  partitioning  of  the  data  set  using 
EM  algorithm  and  the  function  PAPEX  to  determine  the  transformation  of  the  orthogonal  set  of  eigenvalues 
and  eigenvectors  of  the  auto-covariance  matrix  among  the  networks.  The  output  of  an  integrating  unit  in  our 
committee  machine  is  a  transformational  matrix  of  the  image  pair. 

Based  on  those  correspondences  established  by  applying  the  transformational  matrix  between  the  two  images, 
the  control  points  of  two  images  are  obtained.  The  final  step  is  to  recover  the  deformation  of  the  object  by 
calculating  the  polynomial  transformation  using  piece-wise  interpolation.  The  goal  is  to  fine  tune  the  alignment 


927 


Proc.  SP1E  Vol.  4684 


r 


achieved  in  the  initial  registration  by  considering  the  object  as  a  non-rigid  body.  This  allows  for  the  consideration 
of  the  non-rigid  deformation  between  the  images.  The  polynomial  based  transform  such  as  thin-plate  spline 
(TPS)  has  been  shown  to  be  able  to  handle  non-rigid  deformations  of  the  object.  In  our  case,  we  develop  a 
three-layer  MLP  neural  network  to  acquire  the  polynomial  transform,  where  the  input  nodes  axe  fed  with  the 
coordinates  of  the  corresponding  control  points  of  the  image.  The  MLP  is  trained  by  the  corresponding  control 
point  pairs  as  its  inputs  and  outputs  using  the  backpropagation  algorithm.  The  three-layer  MLP  converges  to 
a  polynomial  transform  that  captures  the  relationship  between  corresponding  control  point  pairs.  The  MLP 
neural  network  further  recovers  the  underlying  non-linear  deformation  between  two  images  by  interpolating 
the  learned  relationship  between  corresponding  control  point  pairs.  Finally,  by  applying  the  MLP  polynomial 
transform  to  all  the  pixels  in  an  image  followed  by  a  bilinear  grey  level  interpolation  method,  we  can  obtain  all 
the  pixel  intensity  values  in  the  image  that  is  our  final  registered  image.  ' 

In  order  to  apply  our  registration  method,  we  need  to  establish  the  corresponding  objects  between  two 
images.  We  have  developed  a  stochastic  segmentation  method  to  extract  the  corresponding  objects  in  both 
images.  The  mPAR  uses  these  corresponding  objects  to  combine  multiple  transforms  recovered  by  PAR  in  a 
mixture  form.  After  the  mPAR  step,  we  can  easily  establish  control  point  correspondence  where  control  points 
are  selected  comer  points  on  the  boundaries  of  extracted  corresponding  objects.  Next,  we  will  describe  our 
segmentation  method  for  control  point  extraction. 


2.3.  Segmentation 


Assume  that  each  pixel  in  the  image  can  be  decomposed  into  pixel  images  x  and  context  image  1.  By  ignoring 
information  regarding  the  spatial  ordering  of  pixels,  we  can  treat  context  images  as  random  variables  and  describe 
them  using  a  multinomial  distribution  with  unknown  parameters  7r*,  k  =  1, ...,  K ,  which  can  be  interpreted  as 
a  prior  probability  of  pixel  labels  determined  by  the  global  context  information.  In  particular,  based  on  the 
statistical  properties  of  pixel  images,  where  pixel  image  is  defined  as  the  observed  gray  level  associated  with  the 
pixel,  a  SFNM  distribution  is  justified  to  model  the  image  histogram  by  determining  the  optimal  parameters 
with  respect  to  a  distance  measure  of  a  sum  of  the  following  general  form  /r(z)  —  SfcLi  'jrk9(z\fJ‘k>crk)  with 

*  #  » j 

£  7Tk  =  1,  Jr*  >  0,  and  g(x\/xk:al)  =  -jrrr  where  / ik  and  ak  axe  the  mean  and  variance  of 


k=l 


the  kth  Gaussian  kernel,  and  7vk  is  the  global  regularization  parameter.  We  use  K  to  denote  the  number  of 
Gaussian  components  and  r  to  denote  the  parameter  vector.  This  tissue  quantification  is  achieved  through 
three  completely  unsupervised  steps:  (1)  parameter  initialization  by  optimal  histogram  quantization,  (2)  model 
estimation  by  histogram  based  fast  EM  algorithm,  and  (3)  model  selection  by  minimum  description  length 
(MDL)  criterion.  It  is  shown  that  the  SFNM  model  converges  to  the  true  distribution  when  the  pixel  image  are 
asymptotically  independent.17 


After  we  obtain  the  optimal  parameters  for  all  components,  we  use  a  multiple  thresholding  procedure  for 
initializing  image  segmentation  based  on  maximum  likelihood  (ML)  principle,  which  is  followed  by  contextual 
Bayes  relaxation  labeling  (CBRL)  algorithm  to  obtain  a  consistent  labeling  solution  based  on  localized  SFNM 
formulation  for  improving  initial  segmentation  by  using  neighborhood  text  regularities.  We  define  the  component 
in  localized  SFNM  by  the  support  function: 


SiW=’‘)^“p(-i£2^r)  (22) 

where  7 rj^  is  the  local  conditional  prior  of  regions,  the  support  function  Si(k)  is  a  function  of  the  component 
(tissue  type)  k.  Then  tissue  segmentation  is  interpreted  as  the  satisfaction  of  a  system  of  inequalities  as  follows 
>  Si(k)  for  all  k  and  for  i  =  1, ...,  IV,  where  =  1, ...,  N  is  the  context  images,  and  a  consistent  labeling 
is  defined  as  the  one  having  maximum  support  at  each  pixel  simultaneously.  We  further  define  the  average 

N 

local  consistency  measure  A(l)  =  I(li,  k)Si(k)  to  link  consistent  labeling  to  global  optimization.18 

i=l 
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Figure  1:  Segmentation  result  of  fibroglandular  tissue  (left:  pre-contrast  image;  right:  post  contrast  image). 


Figure  2:  Extracted  control  points  in  a  pre-contrast  image  (left)  and  a  post-contrast  image  (right). 

3.  RESULT  AND  DISCUSSION 

We  applied  out  mPAR  method  to  register  pre-  and  post-contrast  MR  breast  images.  A  3-D  breast  MRI  scan  is 
acquired  prior  to  the  injection  of  contrast  agent,  followed  by  a  sequence  of  3-D  breast  MRI  scan  after  the  agent 
has  been  applied.  The  goal  of  image  registration  for  this  study  is  to  relate  any  point  in  the  post-contrast  sequence 
to  the  pre-contrast  reference  image.  The  first  step  of  our  method  is  try  to  correct  the  global  Tnisaliernmpnt.  by 
performed  a  pre-registration  step  using  the  PAR  method.  The  PAR  method  uses  the  extracted  skinline  as  the 
control  objects  to  globally  align  images.  Since  we  do  not  want  to  change  any  contents  in  the  post-contrast  image, 
we  decide  to  apply  our  registration  method  to  the  pre-contrast  image.  After  we  performed  the  pre-registration 
step,  we  applied  our  segmentation  method  to  extract  the  control  objects  (fibroglandular  tissue)  in  both  pre- 
and  post-contrast  images.  Figure  1  shows  the  segmented  regions  of  fibroglandular  tissue  that  will  be  used  in 
our  mPAR  method  as  the  corresponding  control  objects  for  non-rigid  registration. 

Our  mPAR  method  treats  the  regions  of  fibroglandular  tissue  as  multiple  objects  that  can  be  obtained 
using  the  EM  soft-partitioning  algorithm.  We  then  apply  the  PAPEX  algorithm  to  determine  transformation 
matrices  between  all  the  corresponding  objects.  A  committee  machine  is  then  used  to  combine  individual 
transforms  to  form  a  total  transformational  geometry  of  the  non-rigid  object,  i.e.,  resulting  in  a  mixture  of 
individual  transforms.  After  we  apply  the  mixture  transformation  matrix  to  pre-contrast  image,  all  the  objects 
(Le.,  the  regions  of  fibroglandular  tissue)  are  well  aligned.  We  then  apply  a  comer  detection  algorithm  to 
extract  those  high-curvature  points  on  the  boundaries  of  fibroglandular  tissue  to  serve  as  control  points  for  the 
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Figure  3.  A  pair  of  a  contrast-enhanced  breast  MR  images  (left:  pre-contrast  image  with  global  alignment  based  on 
PAR  method;  right:  post  contrast  image). 


Figure  4.  A  pre-contrast  image  (left)  and  a  registered  post-contrast  (right)  image  using  the  mPAR  registration  method. 


MLP.  Figure  2  shown  the  extracted  control  points  in  both  pre-  and  post-contrast  images.  The  control  point 
correspondence  between  pre-  and  post-contrast  images  can  be  conveniently  established  using  a  simple  correlation 
method.  Finally  an  MLP  neural  network  is  trained  using  the  corresponding  control  point  pairs  to  recover  the 
nonlinear  deformation  between  two  images.  The  non-linear  deformation  is  represented  by  a  polynomial  mapping 
function  that  is  interpolated  by  the  neural  network  based  on  the  control  point  correspondence.  Figure  4  shows 
a  registered  pre-  and  post-contrast  image  pair  using  our  non-rigid  registration  method  -  the  mPAR  method. 

To  relate  a  post-contrast  image  to  a  pre-contrast  image,  we  simple  do  the  image  difference  between  the 
post-contrast  image  and  the  pre-contrast  image.  The  difference  results  are  the  enhanced  area  of  fibroglandulax 
tissue,  i.e.,  the  area  that  has  been  highlighted  by  the  contrast  agent.  Figure  5  (left)  shows  the  difference  images 
between  the  post-contrast  image  and  registered  pre-contrast  image.  As  a  comparison,  we  also  show  the  difference 
image  of  the  post-contrast  image  and  unregistered  (original)  pre-contrast  image  in  Figure  5  (right),  where  false 
regions  of  the  enhanced  area  of  fibroglandulax  tissue  can  be  observed  due  to  the  misalignment  between  pre- 
and  post  contrast  images.  The  results  demonstrate  that  the  mPAR  can  successfully  recover  the  deformation 
between  pre-  and  post-contrast  images  so  the  enhanced  area  can  be  accurately  extracted  by  image  differencing 
between  the  post-contrast  image  and  registered  pre-contrast  image. 
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Figure  5.  Difference  images  of  the  enhanced  area  (left:  difference  image  with  our  mPAR  registration;  right:  difference 
image  without  registration). 


4.  CONCLUSION 

In  this  paper,  we  have  presented  the  theoretical  concepts  and  methods  of  a  neural  computation  based  non-rigid 
image  registration  method  -  the  mPAR  method.  The  approach  uses  a  committee  machine  to  recover  the  total 
transformational  geometry  of  the  non-rigid  object.  The  committee  machine  combines  multiple  rigid  transforms 
in  a  finite  mixture  registration  scheme.  Finite  mixture  transform  combination  is  a  novel  technique  to  combine 
multiple  transforms  contained  in  a  single  image,  which  yields  a  lower  mean  square  error  (MSE)  than  that  of  using 
any  local  transform  alone.  In  other  word,  the  mPAR  method  results  in  a  smooth  image  while  local  transforms 
yield  an  image  containing  discontinuities  on  transform  boundaries.  In  addition,  the  registration  obtained  by 
the  committee  machine  is  fine  tuned  using  the  nonlinear  transform  recovered  by  a  MLP  network  The  MLP 
neural  network  is  trained  by  the  extracted  control  point  pairs  and  converges  to  a  polynomial  transform  that 
captures  the  nonlinear  deformation  between  two  images. 

We  have  applied  our  mPAR  algorithm  to  MR  breast  registration  problem.  The  experimental  results  have 
demonstrated  that  our  mPAR  method  can  successfully  recover  the  nonlinear  deformation  between  images. 
More  importantly,  our  mPAR  method  aligns  two  point  sets  without  the  need  of  establishing  explicit  point 
correspondence.  Instead,  it  aligns  the  two  point  sets  by  mfmrmring  the  relative  entropy  between  their  probability 
distributions.  Furthermore,  through  aligning  the  two  point  sets  initially  by  our  mPAR  method,  we  can  improve 
the  accuracy  of  establishing  the  feature  point  correspondence.  The  MLP  neural  network  has  been  introduced 
to  recover  the  deformation  based  on  established  feature  point  correspondence.  The  MLP  has  also  resulted 
in  a  desirable  error-resilience  property  for  control  point  selection  errors.  However,  some  distortion  can  be 
observed  in  the  final  warped  images.  We  believe  that  the  distortion  is  largely  caused  by  the  errors  in  control 
point  selection  and  correspondence.  Improvement  in  this  portion  should  reduce  the  distortion  and  yield  better 
registration  results.  Finally,  it  is  worthy  noting  that  using  neural  networks  in  this  problem  has  increased  the 
generality  of  this  approach  to  allow  us  to  improve  the  algorithm  and  adjust  performance  as  imaging  conditions 
change. 
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